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ABSTRACT 

(SJ " We propose a plausible mechanism to explain the formation of the so-called "obscuring tori" around 

t-H . active galactic nuclei (AGNs) based on three-dimensional hydrodynamic simulations including ra- 

diative feedback from the central source. The X-ray heating and radiation pressure on the gas are 
CN| ' explicitly calculated using a ray-tracing method. This radiation feedback drives a "fountain" , that is, 

i_ rv a vertical circulation of gas in the central a few to tens parsecs. Interaction between the non-steady 

outflows and inflows causes the formation of a geometrically thick torus with internal turbulent mo- 
tion. As a result, the AGN is obscured for a wide range of solid angles. In a quasi-steady state, the 
opening angles for the column density toward a black hole < 10 23 cm -2 are approximately ±30° and 
. ±50° for AGNs with 10% and 1% Eddington luminosity, respectively. Mass inflows through the torus 

£SJ 1 coexist with the outflow and internal turbulent motion, and the average mass accretion rate to the 

central parsec region is 2 x 10~ 4 ~ 10~ 3 A/© yr _1 ; this is about ten times smaller than accretion 
rate required to maintain the AGN luminosity. This implies that relatively luminous AGN activity 
is intrinsically intermittent or that there are other mechanisms, such as stellar energy feedback, that 
enhance the mass accretion to the center. 

Subject headings: galaxies: Seyfert - galaxies: starburst - ISM: structure - ISM: molecules - method: 
numerical 

Of 

2 1 1. INTRODUCTION 

-£3 ■ 

The presence of optically thick obscuring tori has been postulated to explain the various observed properties of active 

■ galactic nuclei (A GNs), particularly with respect to the two major categories of AGNs, namely, type 1 and type 2 
(jAntonucclll993[ ). However, the origin of the optically and geometrically thick material that covers a large solid angle 

. with respect to the AGN and the nature of its three-dimensi onal structure are still u nclear, although its geometry 

■ is often schematically described as a "doughnut-shaped" (e.g.. lUrrv fc Padovani|[r995() . In fact, statistical studies of 
\ the absorption features of local AG Ns suggested have suggested the p resence of complex, inhomogeneous obscuring 

£ — . material around the central engine (|Shi et al.l 120061: iHicks et "ail 120091) . The clumpy torus model is also consistent 
£Nj ■ with the observed spect r al energy distribution s (SEDs) (e.g.. iNenkova et~aT1[200l lElitzunl2"o"06T: iLevenson et"aT1[200l 
\ iSchartmann et al] 120091 : iStalevski et al.l I20121 ). The essential question that then arises is regarding how such an 
inhomogeneous 'torus' is formed and how it attains such a large covering fraction against the energy dissipation of the 
clumpy medium in the gravitational field of the central black hole. 
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Numerous theoretical models have been proposed to explain th e origin of the obs curing material. On an accretion 
disk scale the proposed causes include hydromagnetic disk winds (|Konigl fc Kartidll994D. disk-drive n hydromagnetic 



• • . winds and radiation press ure due to continuum emiss ion from the accretion disk ( Keating et al.ll2012[ ). and warping of 
irradia ted accretion disk s (jArmitage fc Pringldll997l ). Although a pasec-scale dusty torus is observed in some nearby 
k>( ! AGNs (jJaffe et al.ll2004| ). the nucleus could be obscured by relatively less dense ISM in the range of 10-100 pc or even 
by galactic disks themselves (jLevenson et al.lf2001b iGoulding et al.ll2012T ). The re are many instances where evidence 
suggests that most AGNs are associated with circumnuclear starbursts (e.g. Ilmanishi fc Wadal[2004t iDavies et all 
m0% IHicks et al.ll200i IChen et al.l [20091: Ilmanishi et al1[20TTI: IDiamond-Stanic fc Riekd 120121: IWoo et al.ll2012l) . It 
is suggested that the st ar formation activities dri ve mass inflow into AGNs, thereby leading to the growth of black 
holes (see the review by lAlexander fc Hickoxll201lD. The nucle ar sta rbursts themselves could inflate the circumnuclear 
disk and obscure the central source. IWada fc Normanl (|2002l ) and iWada et al.1 (|2009t ) have proposed that a clumpy 
torus-like structure is naturally reproduced due to energy feedback from supernova explosions (SNe) in a disk. In their 
model, the internal density and temperature of the thick disk are highly inhomogeneous, and the velocity field of the 
disk is turbulent. The scale height of the dis k is determined by th e balance between the turbulent energy dissipation 
and energy input from SNe. More recently, Hopkins et al.l (|2012D have proposed that 10-pc-scale thick tori can be 
formed even without the existence of strong stellar energy feedback if gas inflow from kiloparsec scales drives instability 
in the circumnuclear region. 

On a parsec scale or smaller, geometrically and opt ically thick torus sustained by radiation pressure have been sug- 



gested (|krolik fc Begelmanlll988i : iPier fc Kroliklll992D. Fo llowi ng this model , stati c solutions of thick disks supp orted 
by infrared radiation have been explored bv lKrolikl (|2007l ) and lShi fc Krolikl (|2008l ). lOhsuga fc Umemural (|2001l) have 



suggested that static, geometrically thin obscuring walls could be formed by the interplay between the radiation from 
the nucleus and a ring-like nuclear starburst. 
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It would be more natural t o assume that radi ative feedback causes the dynamical evolution of the surrounding ISM 
rather than static structures. iRoth et all (|2012l ) calculated the net rate of momentum deposition due to the radiation 
pressure in t he surrounding gas and e stimated the mass-loss rate by outflow using Monte Carlo radiative transfer 
calculations. iDorodnitsvn et all (|2012[ ) performed '2.5D' (i.e. basic equations are solved in three dimensions with 
axisymmetry) radiative hydrodynamic simulations, and they found that an AGN torus can be better described in terms 
of a ra diationally supported flow rather than a quasi-static obscuring torus. Using 2D simulations, iSchartmann et all 
(|2011[) studied the evolution of dusty clouds irradiated by the AGN and found that the radiative transfer effect is 
significant depending on the column density of the clouds. 

As is usual in multi-dimensional radiation-hydrodynamic simulations, cirtain simplifications have been applied in 
these previous studies, such as assuming flux-limited diffusion approximation and introducing symmetry for the ro- 
tational axis and equatorial plane of the AGN. In this study, we examine the effect of the radiation from the central 
source on the ISM by assuming a different set of simplifications. We performed fully three-dimensional hydrodynamic 
simulations on scale of a few tens pc around an AGN without introducing any symmetries. We took into account 
radiative heating and pressure due to "direct radiation" from the central source along rays for all the 256 3 grid cells; 
however, we ignored the phenomenon of scattering/re-emission. T his is a natural extension of our previous three- 
dimensional hydrodynamic simulation s of tens-of-parsec scale tori (jWada fe Normanl 120021 iWada fe T omisaka 120051 : 
lYamada et al.ll2007HWada et al.ll2009[ ). in which self-gravity of the ISM, radiative cooling, SNe, and UV heating were 
taken into account. Here, we report the first results without SNe feedback. The SED analysis based on the present 
model will appear elsewhere (Schartmann and Wada, in prep.). In this study, we propose a dynamical process leading 
to the formation of thick tori, and we examine their structures in a quasi-steady state. Further, we discuss the changes 
in the column density and the net mass accretion rate depending on the luminosity of the AGN. 



2. NUMERICAL METHODS AND MODELS 
2.1. Numerical Methods 

Based on our three-dimensional, multi-phase hydrodynamic model (|Wada et al.ll2009 t), we included the radiative 
feedback from the central source, i.e., radiation pressure on the dusty gas and the X-ray heating of co ld, warm, and 
hot i onized gas. In contrast to previous analytical and numerica l studies of radiation-dominated tori(e.g.. |Pier fe Krolikl 
IIlMlKrolikl2rMlShi fe Krolikl2008l : lDiamond-Stanic fe Riekdl20Tl IDorodnitsvn et al1[20lTl [20121 ) . we assume neither 
dynamical equilibrium nor geometrical symmetry. Here, we solve fully three-dimensional hydrodynamic equations, 
accounting for r adiative feedback proces ses from the central source using a ray-tracing method. In a manner similar to 
the approach bv lWada fe Normanl (|2002l ). we account for self-gravity of the gas, rad iative cooling for 2 K < T gas < 10 8 
K, uniform UV radiation for photoelectric heating and H2 formation/destruction (jWada et al.l 120091 ): however, we do 



not include the effect of supernova feedback in this study, in order to clarify the effect of the radiation feedback only. 

We plan to investigate the effects of supernova in the circumnuclear starburst in a future study. 

We use a numer ical code based on Eulcrian hydrodynamics with a uniform grid (jWada fe No rman 2001; Wada] [200~i1 : 
IWada et al.ll2009T ). We solve the following equations numerically to simulate the three-dimensional evolution of a 
rotating gas disk in a fixed spherical gravitational potential under the influence of radiation from the central source. 

dp/dt + V- (pv) = 0, (1) 

d{pv)/dt + (v ■ V)v + V P = -p(V$ + f rad ), (2) 

d(pE)/dt + V • [{ P E + p)v] -pvV$> = pr uv (Go) + pT x - p 2 A{T gas , /h 2 , G ), (3) 

VH sg ^4nGp, (4) 

where $(a;) = $ C xtM + ^bhM + &sg{x); p,p,v denote the density, pressure, and velocity of the gas, and the 
specific total energy E = \v\ 2 /2 + p/(j — l)p, with 7 = 5/3. We assume a time- independent external potential 
$ext(r) = -{27/4) 1 / 2 [v 2 /(r 2 + a 2 ) 1 / 2 + v 2 /(r 2 + a 2 ) 1 / 2 ], where a x = 100 pc, a 2 = 2.5 kpc, Vl = 147 km s" 1 , v 2 = 147 
km s"\ and $ B hM = -GM B n/(r 2 + 6 2 ) 1/2 , where M BH = 1-3 x 1O 7 M (see Fig. 1 in IWada et all (pOOl ) for the 
rotation curve). The potential caused by the black hole (BH) is smoothed within r ~ b = 45, where S denotes the 
minimum grid size (= 0.25 pc), in order to avoid too small time steps around the BH. In the central grid cells at 
r < 2(5, physical quantities remains constant. 

We solve the hydrodynamic part of the basic equations using AUSM (Advection Upstream Splitting Method) (|Liou fe Steffenl 
fl993l) . We use 256 3 grid points. The uniform Cartesian grid covers a 64 3 pc 3 region around the galactic center (i.e. 
the spatial resolution is 0.25 pc). The Poisson equation, eq.(|l|) is solved to calculate the self-gravity of the gas using 
the fast Fourier transform and the convolution method with 512 3 grid points along wiht a periodic Gree n's function. 
We s olve the non-equilibrium chemistry of hydrogen molecules along with the hydrodynamic equations ( Wa da et 
[20091) . 

We consider the radial component of the radiation pressure: 

f rad = J ^dv, (5) 
where Xv denotes the total mass extinction coefficient due to dust absorption and Thomson scattering, i.e. \ u = 
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Xdust,u + Xt- The radial component of the flux at the radius r, F r v is 
where t v — J Xvpds. 

The only explicit radiation source used here is an accretion disk whose size is five orders of magnitude smaller than 
the grid size in the present calculations. Therefore we assumed that radiation is emitted from a point source. How ever, 
the radiation flux originating from the accretion disk is not necessarily spherically symmetric (e.g. INetzerlfl987[ ). In 
our study, we simply assume emission from a thin accretion disk without considering the limb darkening effect, i.e., 
L v = 2Lagn \ cos 6* I, where 9 denotes the angle from the rotational axis (z-axis) (see also discussion in §4). 

In the simulations, we considered the radial component of the flux f r rad only for the radiative heating and pressure. 
This approximation is not necessarily correct especially with respect to optically thick dusty gas around the central 
few parsecs; however, beyond this region, the accel eration of the gas d ue to strong X-rays can be expected to be radial. 
Even for the central region spanning a few parsecs. iRoth et all <|2012f ) showed that the gas acceleration is nearly radial 
using three-dimensional Monte Carlo radiative transfer calculations. It is noteworthy that since the hydrodynamics is 
calculated in a fully three-dimensional manner, the radiative acceleration may contribute to gas dynamics along any 
direction. Therefore the re sultant density and velocity structures could be very different from the symmetric ones (cf. 
lOhsuga fe Umemurall200ll) . The optical depth r„ is calculated at every time step along a ray from the central source 
at each grid point, i.e. 256 3 rays are used in the computational box. 

In the energy equation (eq . (3)), we use a cooling function based on a radiative transfer model of photodissociation 
regions (|Meiierink fc Spaansl 120051) . denoted as A(T g , /h 2 , Go), assuming so lar metallicity, wh ich is a function of the 
molecular gas fraction / H 2 and the intensity of far ultraviolet (FUV), G , (Wa da et al.l 12009). The radiative cooling 
rate below ~ 10 4 K is self-consistently modified depending on /h 2 ; which is determined for each grid cell. We include 
the effect of heating due to the photoelectric effect, Tuv- We assume a uniform FUV radiation field, with the Habing 
unit Go = 1000, originating in the star-forming regions in the central several tens parsecs. In reality, both the density 
field and distribution of massive stars should not be uniform, and Go should have a large dispersion with some radial 
dependence. However, the FUV radiation affects the chemistry of the ISM (mostl y distribution of m olecular clouds) 
although this does not significantly alter the dynamics of the circumnuclear ISM (I Wada et al.l 120091) . The effects of 
the non-uniform radiation field with multiple sources on the gas dynamics is an important problem, but it is beyond 
the scope of this paper. 

The temperature of the interstell ar dust, XH„ s t,, at a g iven position irradiated by the central source is calculated by 
assuming thermal equilibrium (e.e. iNenkova et aDl2008f ). and if Td us t > 1500 K, we assume that the dust is sublimate. 
Here we assume that the ISM in the central tens parsecs is optically thin with resp ect to the re-em i ssion from the 
host dust. The SED of the AGN and the dust absorption cross sect ion are taken fromlLaor fc Drainel (|1993l ). We use 
a standard galactic dust model, which is the same as that used bvSchart mann et al. (2010), with eight logarithmic 
frequency bins between 0.001 fim and 100 /im. 

The effects hard X-ray on the therm al and chemical structures of the ISM are significant in the region spanning 
subparsec distances to tens of parsecs fPcrcz-B eaupuits et al.l [201 lh . Further, X-ray heating may be essential for 
galactic scale feedback (jHambrick et al.ll2011[ ). Here we consider heating by the X-ray, i.e. interactions between the 
high ene rgy secondary electrons and thermal electrons (Coulomb heating ), and H2 ionizing heating in the warm and 
cold gas (|Malonev. Hollenbach. fc Tielensll 19961 iMeiierink fc Spa ans 2005^ . The Coulomb heating rate Tx.c for the gas 
with the number density n is given by 

Tx.c = V nH x, (7) 
where r\ denotes the heating efficiency (Dalgarno et al.lll999HMeiierink fc Spaans!l200"5l ). and the X -ray energy deposi- 
tion rate iJx = 3.8 x 10~ 25 £ ergs s _1 and the ionizing parameter £ = AirFx/n — L^e" J T " du /nr 2 , with Lx — O.ILagn- 
For optically thin hot gas with T > 10 4 K, the eff ects of Compton heating and X-ray photoionization heating are 
considered, and the net heating rate (|Blondinlll994[ ) is approximately 

r x , h = 8.9 x 10~ 36 £(Tx - 4T) + 1.5 x 10~ 21 £ 1/4 T- 1/2 (1 - T/T x ) ergs" 1 cm 3 , (8) 
with the characteristic temperature of the X-ray radiation, Tx — 10 8 K. 

2.2. Initial conditions and model parameters 

In order to prepare a quasi-steady initial condition without radiation feedback, we first run a model of an axisym- 
metric and rotationally supported thin disk with a uniform density profile (thickness of 2.5 pc) and a total gas mass 
of M g — 6.68 x 10 6 M©. Random density fluctuations, which are less than 1% of the unperturbed values, are added to 
the initial gas disk. After t ~ 1.7 Myr, the thin disk settles in a flared shape with the total mass 6.64 x 10 6 M Q . Since 
we allow outflows from the computational boundaries, the total gas mass decreases during the evolution. 

A free parameter in the present simulations is the luminosity of the AGN. Here, we explore models with two different 
Eddington ratios Lagn I Le = 0.1 and 0.01, where the Eddington luminosity Le = 47tGcto )9 A/bh/o't. In each model, 
^agn stays constant during calculations, i.e. ~ 3 MyrtQ. 

1 The numerical calculations were carried out on the NEC SX-9 (16 vector processors/1 TB memory/1.6 TFlops) at the National 
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3. RESULTS 
3.1. Onset of Thick Torus 

Figure [T] shows the onset of biconical outflows, cavities, and thick disk after the radiation from the central source 
is initiated. Radiation-driven outflows start forming in the inner region, and subsequently propagate outward. At 
t = 1.986 Myr, "back- flows" toward the disk plane appear outside the biconical outflows. These accretion flows 
interact with the disk, thereby leading to the formation of a geometrically thick disk that has a clear boundary with 
p~ 1O 2 M pc~ 3 (Fig. [TJ. 

The interactions of vertical flows with opposite directions are more clearly seen in Fig. [5] (t = 2.146 Myr), which 
shows the distribution of v z . There are several contact surfaces between upward and downward flows. As a result, 
turbulent motions that support the thick geometrical structure of the disk are generated. The complicated distribution 
of flows in the thick disk is also seen in the velocity field at the equatorial plane (left panel of Fig. [2|). 

The downward flows can be naturally expected due to the balance between radiation pressure and the gravity. The 
ratio between the radiation pressure for the dusty gas and the gravity due to the central BH is 

frad = L AGN Xd\cos6\ ^_ Td 

fgrav l AtkGcMbH ' 

where denotes the optical depth for the dusty gas at a given point. Suppose that the gas density is a function of z, 
e.g., n(z) — N co i/^/2Trhe^ z '^ 2h . Fig. [3] shows how the radiation-pressure-dominated region and gravity-dominated 
region are distributed arou nd the AGN for the given gas density. Here, we assume the dust-to-gas mass ratio of 0.03 
(jOhsuga fc Umemural 12001). Lagn — f0 44 ergs -1 , dust extinction coefficient Xd = fO 5 cm 2 g _1 , and the gas density 
at z = as no = 1000 cm~ 3 . The four curves represent the condition f ra d/ fgravi = 1 , with the vertical extent of the 
gas h = 4, 6, 8, and 10 pc, for which radiation pressure dominates in the inner side of the curve, therefore outflows are 
expected. On the other hand, the radiation pressure cannot prevent gas accretion outside the critical lines. In fact, 
accretion flows appear in the gravity-dominated domain as seen in Figs. Q] and [2] Moreover, since the critical lines 
dividing the two domains are not purely radial (Fig. [3]), outflows with large angles from the z-axis may collide with 
the boundary of the cavity. This pushes the boundaries of the outflow outward, and the gases in the outflows change 
their directions and fall towards the center as indicated by the velocity vectors in Fig,[2j This also causes the biconical 
cavities to become non-steady, and the vertical circulation of the flow, or the "fountain" of gas, circulates around the 
AGN on a scale of 10 pc. The kinetic energy of this circulating flow is partially used to generate the random motion 
in the disk and maintain its thickness. In other words, both radiation energy and gravity are the cause of origin of the 
internal turbulent motion in the dislfi 

3.2. Structures and Dynamics in Quasi-steady State 

Figure |4] shows that the density structures in a quasi-steady states differ depending on the luminosity of the central 
source. Here "quasi-steady" means that the global morphology has reached an approximately steady state; that is the 
phase diagram of the gas (Fig. [5]) does not significantly change after this point. The accretion rates to the central 
region are also approximately constant. This is achieved roughly after t = 3.0 — 3.5 Myr (Fig. [9]). Note that since the 
duration of the present simulation is much shorter than the accretion time scale (~ 10 9 Myr), the 'quasi-steady' does 
not necessary mean the long-term stability. 

In the more luminous model (Lagn I^Edd — 0.1, upper panels), inhomogeneous outflows with p ~ 0.1 — 1 Mq pc -3 
and velocity of ~ 100 km s _1 are formed, and diffuse gas (p ~ 0.1 Mq pc -3 ) extends for over 20 pc above the equatorial 
plane. On the other hand, diffuse and relatively stable outflows with p <~ 10~ 2 M© pc~ 3 or less are formed in the less 
luminous model with Lagn / L>Edd = 0.01. The gas temperature distributions of the two models (Fig. [5J show that 
dense spiral arms and clumps in the equatorial plane exhibit a temperature of less than 1000 K, in contrast to the 
"warm" (~ 10 4 ~ 5 K) torus. In the less luminous AGN model, an incresed number of hot gas (> 10 6 K) regions appear 
in the biconical outflow. This is because the radiative cooling is less effective due to reduced gas density in the region, 
and as a consequence, the X-ray is not attenuated. 

Figure |6] shows the effects of the radiative feedback on the phases of the gas. The dominant phase is the one around 
Tgas ~ 8000 K. The dense and cold gas (p > 10 2 Mq pc -3 and T gas < 1000 K) is not affected by the radiative feedback. 
The X-ray heating produces warm (T gas < 1000 K) gas in the outer regions (i.e., high latitude) of the thick disk, as is 
expected from the XDR model ( Pcrcz-B eaupuits et alj|2~011[ ). Shock-heated gas with a temperature of several 10 4 K 
also appears in the same density region. 

The structure of the torus does not show perfect symmetry for the rotational axis and the equatorial plane. Figure 
[7] shows four different vertical slices of density around the rotational axis of the model shown in the upper panel of 
Fig. IU The diffuse part of torus (p ~ 10 2 ) is slightly tilted, thereby indicating that the outlying region of the disk is 
not completely settled. The shell-like outflows also show asymmetric distributions. It is interesting to note that these 
structures are caused by interactions between the radiation from the central symmetric source and the surrounding 

Astronomical Observatory of Japan and Cyberscience Center, Tohoku University. A single typical run for 3 Myr takes approximately 120 
CPU hours. 

2 In a self-gravitating, multi-phase gas disk, gravitational and thermal instabilities themselves can generate and maintain turbulent 
motion (Wada ct al. 2002). In the present model, random motions of the high-density gas near the equatorial plane (z ~ 0) are dominated 
by this gravity-driven turbulence. However, the random motions along the vertical direction in the less dense gas in the thick disk are not 
driven by gravitational instability. 
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Fig. 1. — Initial evolution of radiation-driven outflow and onset of thick disk in a model with ^AGN/^Edd = 0.1. The vertical slices 
represent the x-z planes. Arrows represent gas velocity (a unit vector of 500 km s _1 is shown as reference along the y-axis.) 

inhomogeneous gas, thereby demonstrating importance of the three-dimensional treatment of radiation hydrodynamics 
in understanding the structure of AGN. 
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Fig. 2. — Observed ^-component of velocity at t = 2.15 Myr (Lagn/Le = 0.1). The vertical slice indicates the x-z plane. The downward 
flows (i.e., flow heading to the equatorial plane) appear outside the region of the biconical outflow, and these flows interact with the opposite 
directed flows, which are represented by the regions indicated in green/blue and red/yellow. The upward and downward flows coexist and 
generate turbulent motions in the equatorial plane (left panel) as well as in the vertical direction. 




x [pc] 

Fig. 3. — Critical lines between radiation-pressure dominant and gravity dominant regions for gas disk with different vertical extentions, 
represented by n(z) = no/(v / 2~7r/i)e -z ^ 2h with h = 4, 6, 8, and 10 pc. The dashed line indicates the critical line for h = 4 pc and 
Lagn = 10 46 erg s — 1 , which is 100 times that for the other cases. 

Obscuration due to the vertically extended gas in the central region over tens of parsecs and how it differs between 
models are shown in Fig. [51 which shows the plot of the column density toward the central BH as a function of the 
viewing angle. If there is no radiation from the central source, the opening angle for N < 10 23 cm~ 2 is approximately 
140°. The radiation feedback causes the thickness of the disk to increse, and the opening angles are approxiamtely 60° 
and 100° for Lagn/Le = 0.1 and 0.01, respectively. It is noteworthly that the column density does not only depend 
on the viewing angle, but also on the azimuthal angles. The column density varies by a factor of ten or more for a 
given viewing angle, thereby indicating the internal inhomogeneous structure of the disk. The structure of the disk is 
similar to that generated by supernova explosions (jWada fc Normanll2002t IWada et al.l [20091. However, the present 
simulations imply that the geometrically and optically thick torus can be naturally formed due to gravitational energy 
in the central part with the aid of the radiation feedback in the case of the AGNs with low Eddington ratios. 

The radiation pressure caused by the central source usually has "negative" effect on the mass accretion toward the 
center; however, our results indicate that mass accretion is not terminated by radiative feedback. Figure [9] shows the 
time evolution of the mass of the two models in four radial bins. The lower luminosity model [Lagn I^Edd = 0.01 ) 
shows continuous net mass accretion in the central 8-pc region, and there is almost no net inflow beyond r ~ 4 pc. 
The mass accretion in the central part is caused by gravitational torque due to the non-axisymmetric features near 
the equatorial plane of the central part of the disk (Fig. [I]). 

3 The amount of the angular momentum transfer during one rotational period by m-armed spiral density waves is roughly proportional 
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Fig. 4. — Gas density in quasi-steady state of two models: (top) Lagn/Le = 0.1 at t = 4.55 Myr and (bottom) Lagn/Le = 0.01 at 
t = 4.59 Myr. The vertical slices indicate the x-z planes. . . ... _ o _ q 

Ihc accretion rate to the central parsec increases with decreasing AGN luminosity: Al = 1.7 x 10 M Q yr and 
M = 1.3 x lO~ 4 M yr^ 1 for L AGN /L E = 0.01 and 0.1, respectively. For Lagn = 0.1L B , the net accretion rate is 
one order of magnitude smaller than the accretion rate required to constant AGN luminosity; the AGN luminosity is 
assumed to be constant during the calculations provided that the energy conversion efficiency is 0.1. On the other hand, 
the observed accretion rate in the low luminosity AGN model is comparable to that required for the luminosity of the 
source. These results may suggest that luminous AGNs are intermittent, and their luminosity is not instantaneously 
determined by mass inflow from r ~ 1 pc. 

4. CONCLUSION AND DISCUSSION 

In our study, we showed that a geometrically and optically thick torus can be naturally formed in the central region 
extending tens parsecs around a low luminosity AGN. We found that radiation drives a "fountain" of gas, i.e. vertical 
circulation, and this also naturally produces a thick, turbulent, and inhomogeneous disk that resembles a torus. The 
opening angle of the torus for N < 10 23 cm~ 3 is 60 — 100° for Lagn = 0.1 — O.OILe, in contrast to 140° for the model 
without radiative feedback. The average accretion rate to the central parsec region is larger in a less luminous AGN; 
however, this accretion rate is insufficient to maintain the AGN luminosity over several Myrs. These results imply that 
the AGN luminosity and structure of the surrounding ISM alternate between an active phase (high luminosity and a 
torus with a small opening angle) and inactive phase (low luminosity and a torus with a large opening angle). For 
thick tori (Fig. U) to maintained their structures for several million years or longer, there should be other mec hanisms 
to enhance the accretion rate, such as supernova-driven turbulence (|Wada fc Normanll2002| : iWada et aLll2009f ). stellar 
mass loss (ISchartmann et al. I2010D . turbulence drives by galactic inflows (Hopkins et al.l [2012D andradiation drag 



(jKawakatu fc Umemurall2002 ) . In fact, recent observatio ns suggest a tight correlation between star formation rate and 
the Eddington ratio of the central BH (jChen et al.ll2009f ). 

Figure [8] shows that the covering fraction of the central source is smaller in AGNs with lower luminosity. This is 
because more mass can be ejected from the central region to higher latitudes, and as a result, the radiation from the 
AGN is self-obscured. Therefore, the radiation-driven outflows appear only in a small solid angle around the rotational 
axis, as explained in §3.1 and illustrated in Fig. [3j This means that more luminous AGNs are obscured over larger soli d 
angles, altho ugh this may ap pear inconsistent with the observed "receding torus" (|Lawrencelll991t lUeda et aLi r2003). 
For example. iHasingerl (|2008| ) studied 1290 AGNs and found that the absorbed fraction decreases with X-ray luminosity 
in the range of Lx = 10 42-46 erg s" 1 . However, it is to be noted that the AGNs explored in the present simulation 

to n ta,n 9 iA^/rn (e.g. Binncy & Tremainc 2008), where A m is a ratio of the amplitude of the spiral to the average density, and 9i is the 
pitch angle of the spiral. For non-linearly developed spiral waves (i.e. A m > 1) with tanfi^ ~ 0.1, m ~ 1, about 3% of the initial angular 
momentum of the disk is transferred in a rotational period. 
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Fig. 5. — Gas temperature in quasi-steady state of two models: (top) Lagn/Le = 0.1 at t = 4.55 Myr and (bottom) Lagn/Le = 0.01 
at t = 4.59 Myr. The vertical slices indicate the x-z planes. 

have relatively low luminosity, i.e., in the range of Lx = 10 42 ~ 43 erg s _1 . In this small luminosity range, the structure 
of the torus does not vary significantly, especially for higher values of column density (N > 10 24 cm -2 ). For a central 
source that is 100 times brighter, the radiation-dominant region is considerably larger, as shown by th e dashed line in 
Fig. [31 and therefore, we can expect a larger opening angle for this case. More recently. ILu et all (|2010D have reportred 
results that are significantly different from those of previous studies; in their SDSS/FIRST survey, they found that the 
type-1 fraction exhibits a constant of 20% in the [O III] 5007 luminosity range of 40.7 < log(L[o ln ]/ergs _1 ) < 43.5. 
They also found that the type-1 fraction is independent of the Eddington ratio for its value between 0.01 and 1, if only 
high density (> 10M Q pc -3 ) gas in the torus contributes to the obscuration; this is also the case as seen in Fig. |4j 

Further, this study shows that the global structures of the torus and outflows are not static, and there is never 
a perfect symmetry in terms of the equatorial plane and the rotation axis (Fig. [7]), even if we assume the presence 
of angle-dependent, axisymmetric flux from the central source. The opening angle of the torus and other structures 
could be affected by the choice of the central radiation source. In this study, we assumed a thin accretion disk; 
however, recent two-dimensional radiation-magnetohydrodynamic simulations have shown that geometrically thick 
and radiatively inefficient a ccretion flows associated with outflows are formed depending on the central gas density 
(jOhsuga fc Mineshigd I2IT1TI ) . In this case, the radiation flux is collimated around the rotation axis, and the effect 
of the radiation on the disk might be less than that observed in the models used in this study. Even for a thin 
accretion disk, the rotational axis is not necessarily aligned to the galactic rotational axis. This causes anisotropic 
illumina tion of the ISM, and this misalignment could also affect the occultation of the central source by the torus 
(see also iKawaguchi fc Morill2011[) . Consequently, we can expect the formation of more such anisotropic structures of 
the torus and outflows in such cases. This can be confir med by high-resol ution observations of nuclear active galaxies 
associated with molecular outflows, such as NGC 1377 (jAalto et al.ll2012T ) and by ALMA in the near future. 
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NAOJ and Cyberscience Center, Tohoku University. The author thanks N. Horiuchi, Y. Sakaguchi, and the NEC 
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Fig. 8. — Total column density distribution as function of viewing angles (90°: polc-on view and 0°: edge-on view) and azimuthal angles. 
For each viewing angle, the column density is plotted for 50 azimuthal angles. The solid line indicates the azimuthal averaged column 
density for the viewing angles, (top) Before the radiation feedback is turned on (t = 1.67 Myr), (middle) LagnI^e = 0-01 at t = 4.55 
Myr, (bottom) L AGN /L E = 0.1 at t = 4.68 Myr. 
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Fig. 9. — Evolution of radial mass, (left) model qi (Eddington ratio = 0.1) and (right) model qg (Eddington ratio = 0.01) for r < 1.0 pc 
(thick line), 1 < r < 2 pc (dotted line), 2 < r < 4 pc (dashed line), and 4 < r < 8 pc (dot-dashed line). 



